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ABSTRACT 

, We investigate the spins and shapes of over a million dark matter haloes identified 

^ • at z = in the Millennium simulation. Our sample spans halo masses ranging from 

f — | dwarf galaxies to rich galaxy clusters. The very large dynamic range of this ACDM cos- 

. mological simulation enables the distribution of spins and shapes and their variation 

^sD ■ with halo mass and environment to be characterised with unprecedented precision. We 

compare results for haloes identified using three different algorithms, and investigate 
^— ^ ■ (and remove) biases in the estimate of angular momentum introduced by both the al- 

gorithm itself and by numerical effects. We introduce a novel halo definition called the 
TREE halo, based on the branches of the halo merger trees, which is more appropriate 
• for comparison with real astronomical objects than the traditional "friends-of-friends" 

Qh' and "spherical overdensity" algorithms. We find that for this many objects, the tra- 

q . ditional lognormal function is no longer an adequate description of the distribution, 

' P(A), of the dimensionlcss spin parameter A, and we provide a different function 

! that gives a better fit for TREE and spherical overdensity haloes. The variation of 

C$ • spin with halo mass is weak but detectable, although the trend depends strongly on 

' the halo definition used. For the entire population of haloes, we find median values 

.,— I [ of Amed = 0.0367-0.0429, depending on the definition of a halo. The haloes exhibit 

■ a range of shapes, with a preference for prolateness over oblateness. More massive 

\ haloes tend to be less spherical and more prolate. We find that the more spherical 

haloes have less coherent rotation in the median, and those closest to spherical have a 
spin independent of mass (A me d ~ 0.033). The most massive have a spin independent 
of shape (A me d ~ 0.032). The majority of haloes have their angular momentum vector 
aligned with their minor axis and perpendicular to their major axis. We find a gen- 
eral trend for higher spin haloes to be more clustered, with a stronger effect for more 
massive haloes. For galaxy cluster haloes, this can be larger than a factor of ~ 2. 

Key words: cosmology: dark matter - galaxies: haloes - methods: TV-body simula- 
tions 
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1 INTRODUCTION 

The formation of galaxies is intimately linked to the acqui- 
sition and distribution of angular momentum. In the cur- 
rent cosmological paradigm, the inflationary ACDM model, 
cosmic structures grow hierarchically. Dark matter haloes 
form by the dissipationless gravitational collapse of mate- 
rial associated with peaks in the primordial density fluctua- 
tion field, growing, through mergers and smooth accretion, 
into objects with a wide range of masses at the present day. 
Galaxies form when baryons cool and condense near the 
centre of these haloes (White & Rees 1978, White & Frenk 



1991). They undergo mergers and tidal interactions along 
with their haloes, giving rise to the rich spectrum of galaxy 
types and environments that we see today. 

Understanding the generation and evolution of the an- 
gular momentum of dark matter haloes is a prerequisite for 
understanding the angular momentum and morphology of 
galaxies. For example, the distribution of halo spins is a 
basic input to models of galaxy formation (e.g. Mo et al. 
1998, van den Bosch 1998, Cole et al. 2000). The early evo- 
lution of the angular momentum of a density perturbation is 
adequately described by the linear tidal torque theory (see 
e.g. Hoyle 1949, Peebles 1969, Doroshkevich 1970, White 
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1984, Catelan & Theuns 1996, Lee & Pen 2000). However, as 
work such as that of White (1984) and Porciani et al. (2002) 
has shown, the non-linear effects inherent in the formation 
of large-scale structure lead to large quantitative disagree- 
ments between the predictions of the tidal torque theory and 
the angular momenta found in iV-body simulations of dark 
matter haloes. 

iV-body simulations provide the way to progress beyond 
the linear regime. As computing power has improved, so has 
the scale and resolution of simulations. Very early numerical 
studies of the angular momentum of "proto-galaxies" were 
performed by Peebles (1971) (with JV P ~ 100 particles) and 
Efstathiou & Jones (1979) (N p = 1000), and led the way to 
the analysis of the spins and shapes of CDM haloes in more 
sophisticated simulations (Davis et al. 1985; Barnes & Efs- 
tathiou 1987; Frenk et al. 1988; all with JV P = 32 768). War- 
ren et al. (1992) used a much larger simulation (jV p ~ 10 6 ) 
and focused particularly on the details of the distributions 
of halo spins and shapes, and their relationship through the 
alignment of the halo angular momentum vector. Cole & 
Lacey (1996) also investigated the shapes and spins of dark 
matter haloes, in addition to various other aspects of halo 
structure. 

An early study by Barnes & Efstathiou (1987) examined 
the relationship between spin and the spatial clustering of 
haloes, as measured by the two-point correlation function. 
Later, Lemson & Kauffmann (1999) examined the environ- 
mental dependence of halo properties, and found no cor- 
relations with halo spin. More recently, Faltenbacher et al. 
(2002) have carried out a mark correlation function analy- 
sis to investigate how spin varies with halo pair-separation. 
They found that neighbouring cluster pairs tend to have 
higher spins than the average. 

Recent years have seen a large amount of work on 
the analysis of haloes in ACDM simulations. Halo shapes 
and their variation with mass were investigated by Bullock 
(2002), Kasun & Evrard (2005), Shaw et al. (2006) and All- 
good et al. (2006). In agreement with previous studies, halo 
spin was found to vary little, if at all, with halo mass. The re- 
lationship between halo shape and spin was investigated by 
Bailin & Steinmetz (2005), Avila-Reese et al. (2005), Shaw 
et al. (2006) and Gottlober & Turchaninov (2006). While 
this paper was being completed, independent analysis of 
halo properties, investigating halo concentrations, spins and 
shapes in a series of simulations, was posted by Maccio et al. 
(2006). 

Using the the 10-billion particle Millennium simulation 
of the evolution of dark matter in the ACDM cosmology 
(Springel et al. 2005), we re-examine, in this paper, some of 
the shape and spin properties of dark matter haloes previ- 
ously considered. Our analysis improves upon earlier work 
because the millions of haloes that formed in this simulation 
provide unprecedented statistical power. This allows us, for 
example, to quantify the distribution of halo spins and the 
relationship between spin, halo mass and shape with a pre- 
cision that has not hitherto been possible. Unlike previous 
work, we consider different ways to identify haloes in the 
simulations; it turns out that the details of halo definition 
and selection can have a strong impact on the results. Fi- 
nally we investigate how halo clustering depends on spin and 
shape for haloes of different masses. 

This paper is structured as follows. Section 2 provides 
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Table 1 . Cosmological and simulation parameters for the Millen- 
nium Run. The first row describes the cosmology used. It gives 
the density parameters Q; := p Pi , where the critical density 

<° crit := IS? ( with P A := 8^G and Otot = + H M = 1), the 
Hubble parameter where Ho = 100/i(km s" 1 Mpc - 1 ), the spec- 
tral index, n, and ag, the linear-theory power spectrum variance 
in spheres of radius 8/i _1 Mpc. The second row gives the simu- 
lation parameters: the length of the side of the simulation cube, 
the number of particles, the resulting particle mass and the grav- 
itational softening length. For further details, see Springel et al. 
(2005). 



a description of the Millennium simulation itself, and the 
various halo properties we shall be investigating. Section 3 
describes the construction of the catalogues whose haloes we 
investigate, including the group-finding algorithms and halo 
selection criteria. These we use to remove haloes whose prop- 
erties are unreliable or biased, due to both numerical effects 
and the group-finding algorithms themselves. The main re- 
sults of this paper are presented in §4, where we describe the 
distribution of halo spins as a function of mass and shape, 
and examine its effect on halo clustering. Our conclusions 
are presented in §5. Finally, the Appendix shows various 
examples of haloes that illustrate the effects of the group- 
finders. 



2 DARK MATTER HALO PROPERTIES IN 
THE MILLENNIUM SIMULATION 

2.1 The simulation 

The Millennium simulation is described by Springel et al. 
(2005). It followed the evolution of 10 billion dark mat- 
ter particles in the ACDM model, the standard paradigm 
of modern cosmology. This is is strongly favoured by mea- 
surements of the temperature anisotropies in the microwave 
background radiation (Spergelet al. 2003, 2006) and by mea- 
surements of the clustering of galaxies (Percival et al. 2002; 
Tegmark et al. 2004; Sanchez et al. 2006). For reference, the 
key parameters of the simulation are listed in Table 1. 

The simulation gives the positions and velocities of all 
particles at each "snapshot" in time; 64 snapshots are stored 
from the initial redshift (z = 127) to the present, enabling 
redshift-dependent statistics to be studied. 

2.2 Halo properties 

For each halo-finding algorithm we consider in this paper, we 
compute a range of halo properties. We shall discuss these 
quantities here, and defer the details of the halo-finding and 
halo-selection algorithms to the next section. 

Much of this work concentrates on the dimensionless 
spin parameter A, introduced by Peebles (1969, 1971). This 
is defined as: 
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where Mh is the halo mass, J is the magnitude of the angular 
momentum vector J (and j is the specific angular momen- 
tum), E is the total energy, and G is Newton's gravitational 
constant. It is important to note that A is defined for any 
object which has a well-defined j, E and Mh; all these quan- 
tities are conserved for an isolated system, virialised or not. 
In reality, none of the haloes in cosmological simulations are 
completely isolated, leading to ambiguities in the definition. 
This means that for the A of a halo to be useful, it is the def- 
inition of halo that requires the most care (and conditions 
such as virialisation; see section 3.2), not the definition of A. 

The meaning of A is therefore best understood by con- 
sidering an isolated, virialised, spherical system. The spin 
parameter can be seen to be a measure of the amount of 
coherent rotation in a system compared to random motions. 
For a spherical object, it is approximately the ratio of its 
own angular velocity to the angular velocity needed for it 
to be supported against gravity solely by rotation (see e.g. 
Padmanabhan 1993). 

The specific angular momentum j and kinetic energy T 
of each halo containing N p particles are given by: 
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where n is the position vector of particle i relative to the 
halo centre, and Vi is its velocity relative to the halo centre 
of momentum. 

The halo potential energy, U, is calculated using all halo 
particles if Mh ^ 1000m p , and is rescaled up from that of 
1000 randomly-sampled particles otherwise. The potential is 
that used in the simulation itself: 
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where JV sel is the number of selected particles (iV sc i ^ 1000), 
r) is the softening length (see Table 1), ry is the magnitude 
of the separation vector between the ith and jth particles in 
the halo, and the softening kernel (see Springel et al. 2001) 



equivalent to rotating the coordinate system to find a set of 
axes in which a torque about one does not induce a rotation 
about another; i.e. such that J is parallel to uj. These axes 
then describe a hypothetical uniform ellipsoid whose axes 
a ^ b ^ c are those of the halo itself: 
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The eigenvalues are the moments of inertia I for rotation 
about that axis. For example, rotation about the semimajor 
axis a has the moment of inertia T a = |Mh(& 2 + c 2 ); note 
that l a ^Z c - These eigenvalues can then be combined 
to find the relative axis lengths, e.g. 
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The axis vectors are given directly by the corresponding 
eigenvectors, so that, for example, rotation about the a-axis 
of the ellipsoid (with the longest length a) has the smallest 
moment of inertia, T a - 

Much of the literature on halo shapes uses the following 
description of the mass distribution, confusingly also calling 
it the inertia tensor (see, e.g. Cole & Lacey 1996, Hopkins 
et al. 2005, Bailin & Steinmetz 2005, Shaw et al. 2006): 
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The results are entirely equivalent: if one diagonalises this 
matrix, then a, b and c can be found as just the square roots 
of the eigenvalues, and the eigenvectors again give the axis 
vectors. 

Once the halo's principal axes have been found, rela- 
tionships between the axes and between the shape and other 
properties such as spin can be examined in terms of the 
axis ratios p := c/b, q := b/a and s := c/a. The minor-to- 
major axis ratio s is a useful measure of the sphericity of 
the system, but does not specify in what way a halo might 
be aspherical. For this, we can use the triaxiality parameter 
introduced by Franx et al. (1991): 
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This measures whether a halo is prolate (T = 1) or oblate 
(T = 0), but it does not quantify how aspherical a halo is. 
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A halo's shape is derived from its mass distribution, 
which we characterise using the inertia tensor, I. This relates 
angular momentum J and angular velocity w through J = 

and is formed from the following components: 



la/3 ■— mi 



(6) 



where ri is the position vector of the ith particle, a and 
(3 are are the tensor indices with values of 1, 2 or 3, and 
5 a is the Kronecker delta. The process of diagonalising I is 



3 THE HALO CATALOGUES 
3.1 Groupfinder algorithms 

The problem of how best to identify groups of particles 
within iV-body simulations is a ubiquitous feature in stud- 
ies of dark-matter halo properties. Many solutions have 
been found to this problem, typically involving first find- 
ing candidate-halo centres, followed by an iterative scheme 
to shrink or grow the halo according to criteria involv- 
ing the binding energy or overdensity. Such algorithms in- 
clude Bound Density Maxima (BDM, Klypin et al. 1999), 
Spline Kernel Interpolative DENMAX (SKID, Weinberg 
et al. 1997), the AMIGA Halo Finder (AHF, Gill et al. 2004), 
and Subfind (Springel et al. 2001), which we use below. 
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The very simple 'friends-of-friends' group-finder (FOF, 
Davis et al. 1985) was run on-the-fly, during the Mil- 
lennium simulation run, with a linking length of so = 
K^box/^part) 1 '' 3 where b = 0.2, to attempt to select viri- 
alised structures in the particle distribution. As is often the 
case (e.g. Porciani et al. 2002, Shaw et al. 2006, Maccio 
et al. 2006), this simple FOF catalogue forms the basis for 
the more sophisticated halo definitions we use. 

An enhanced version of the Subfind program (Springel 
et al. 2005) was run on the data to identify self-bound sub- 
structures within each FOF halo, which we then use to con- 
struct the different group catalogues we investigate. The 
Subfind algorithm is essentially a two-step process. The 
first task is to identify subhalo candidates within each FOF 
halo. This is done using an adaptively-smoothed dark mat- 
ter density field, effectively lowering a density threshold and 
identifying the peaks that grow out of it. The second stage 
consists of performing an iterative gravitational unbinding 
procedure on the candidates, successively removing parti- 
cles that are not bound to the subhalo candidate. For this 
purpose, the potential energies are computed using a tree 
algorithm similar to that used for the simulation itself. The 
candidates that are left with at least 20 particles after this 
procedure are then subhaloes of the parent halo. The algo- 
rithm can and does identify subhaloes within subhaloes. It 
results in the FOF haloes typically consisting of a hierarchy 
of self-bound structures (which are not necessarily bound to 
each other), and a set of particles (referred to as "fuzz") that 
are spatially linked to the halo but not part of any self-bound 
(sub) structure. The most massive "subhalo" (MMSH) typ- 
ically contains most of the mass of the corresponding FOF 
object, and so is best regarded as the self-bound background 
halo itself, with the remaining subhaloes as its substructure. 

In addition to finding the bound structures within 
haloes, Subfind also computes certain subhalo properties, 
which are then stored in the subhalo catalogue files. These 
include the location of the potential minimum, the ID num- 
ber of the most bound particle, the mass (number of parti- 
cles), and the half-mass radius. Subfind also computes and 
stores additional data related to each parent halo. Starting 
at the potential minimum of the MMSH of a halo, three radii 
are found: the first two are those where the density within 
them drops below 200p cr it and 200p, to aid comparison with 
other work. The third is the virial radius proper, which uses 
the fitting formula of Bryan & Norman (1998) for spherical 
top-hat collapse in a flat (Q\ + Qm = 1) cosmology (see also 
Eke et al. 1996): 

— = 18tt 2 + 82(f2 M (z) - 1) - 39(fi M (z) - l) 2 (11) 

Pcrit 

This gives p/p cr it ~ 94 at z — 0. Although these properties 
are associated with sets of grouped particles, Subfind does 
not restrict itself to these particles when computing them. 

We will now describe the group catalogues whose halo 
properties we have investigated, and how they are built from 
the results of the FOF and Subfind algorithms. A key point 
which is used for each halo definition is that we take the halo 
centre to be at the potential minimum of the MMSH. For 
reference, Table 2 gives a list of each halo catalogue we will 
discuss, and the number of haloes they contain. 



Hello Catalogue 


Number of haloes 


FOF (without Subfind) 


17 709121 


('raw' catalogue) 


FOFall 


15 494 624 


(87.5% of FOF raw) 


FOFclean 


1 332 239 


(8.60% of FOFall) 


SOall 


15 458 379 


(99.8% of FOFall) 


SOclean 


1239 494 


(8.02% of SOall) 


TREEall 


17041498 


(110.0% of FOFall) 


TREEclean 


1 503 922 


(8.83% of TREEall) 



Table 2. Numbers of haloes in the halo catalogues defined in 
this paper. The "all" catalogues are the unfiltercd results of the 
groupfinder algorithms described in sections 3.1.1, 3.1.2 and 3.2.2. 
The "clean" catalogues have been filtered as discussed in sections 
3.2.3 and 3.2.4, with the quasi-cquilibrium parameter Q = 0.5 
and the particle-number limit N p ^ 300. 



3.1.1 FOF haloes 

The basic FOF algorithm, run as already described, re- 
sults in 17 709 121 groups containing at least 20 particles, 
i.e. of mass > 1.7 x 10 10 /i _1 Mq, at redshift z = 0. The 
most massive group contained 4 386 162 particles (w 3.8 x 
10 15 /i -1 Mq). Using this catalogue we can therefore iden- 
tify dark matter haloes over a range of more than 5 orders 
of magnitude in mass, ranging from subgalactic clumps to 
the most massive clusters. 

In practice, since we use the centre (potential minimum) 
of the MMSHs as the centres of the haloes themselves, we 
only use FOF haloes for which Subfind has found bound 
substructures. Haloes without substructure (and hence with- 
out an MMSH) are excluded from our base FOF catalogue. 
This has the effect of reducing the catalogue size by 12.5%, 
preferentially at lower masses. We shall refer to this cata- 
logue as FOFall, and we will not discuss the larger raw FOF 
catalogue further. 



3.1.2 SO haloes 

The properties calculated by Subfind make it possible to 
construct a second halo catalogue, in which each halo con- 
sists of only the particles within i? v j r of the centre of the 
MMSH of the corresponding FOF object (note that these 
particles do not have to be members of the FOF halo). Due 
to the way Subfind constructs the MMSH, this yields haloes 
whose definition is similar to those from a "spherical over- 
density" algorithm (Lacey & Cole 1994), so we will refer 
to them as the SO haloes. We do not impose a lower limit 
on the number of particles comprising these objects; as a 
result, haloes can be identified with masses < 20m p when 
their virial radii encompass fewer particles than their origi- 
nal FOF halo. This is simply a consequence of the algorithm 
employed; later examination of halo spins reveals the need 
for a much higher particle-number limit, as discussed in full 
in Section 3.2.4. This halo definition is similar, but not iden- 
tical, to that used in Maccio et al. (2006). 
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3.2 Better halo catalogues 

3.2.1 Groupfinder problems 

The FOF and SO group-finding algorithms have various 
well-documented drawbacks. Groupfinders such as SO which 
use the overdensity contained within a spherical region tend 
to impose a more spherical geometry on the resulting sys- 
tems. Although this is not a problem for many objects, the 
algorithm can sometimes result in very unnatural-looking 
structures. An example is shown in figure A4. This compares 
a massive FOF halo with the corresponding SO object. The 
centre of mass of the FOF system is well separated from 
the minimum of the potential, and the halo is significantly 
elongated. This results, when growing a sphere around the 
potential minimum to form the SO halo, in the virial over- 
density being reached sooner in one direction than another. 
The SO halo contains particles in low-density regions outside 
the FOF halo, and has a sharp cutoff in another direction 
when the FOF halo continues. The more 'normal' haloes, in 
the background, are much less affected. 

The problems associated with the FOF groupfinder can 
be more extreme, and can affect a greater proportion of the 
haloes. One of the most commonly-cited problems (e.g. Gelb 
& Bertschinger 1994, Governato et al. 1997) is that well- 
resolved objects identified using the FOF algorithm are of- 
ten at risk of becoming linked with neighbouring objects via 
tenuous bridges of particles. Low-mass particle bridges are 
usually extremely transient structures, being just a chance 
grouping of particles at that instant in time. The joining of 
two (or more) otherwise unrelated objects of similar mass in 
this way results in a very large velocity dispersion. Examin- 
ing the halo in velocity-space will clearly show the multiple- 
object nature of the system. An example of a halo formed 
from objects joined by a tenuous bridge can be seen in Fig. 
Al. 

Sometimes more massive haloes can be formed by the 
chaining together of somewhat smaller objects that are un- 
dergoing mergers or close flybys with their neighbours. Their 
multi-object composition can again be seen in velocity space 
as well as in real space, and although their connections are 
likely to be less transitory than in the case of a thin bridge, 
these objects are nevertheless well out of their equilibrium 
state, and so are unhelpful when trying to characterise the 
spins of typical dark matter haloes. See Fig. A2 for an ex- 
ample of a larger multi-object halo. 

A similar effect is that of velocity contamination of 
small objects due to their proximity to more massive ones 
(see Fig. A3 for an example). Just as particles can form 
a bridge between passing haloes at the moment of the 
snapshot, so an individual particle orbit can take it within 
the linking length of a neighbouring halo, without forming 
enough of a bridge for the haloes themselves to be joined. 
The smaller halo will be contaminated by these interloper 
particles, which will have a quite different mean velocity to 
the halo's own particles. This causes the mean velocity to 
be shifted away from that of the 'original' halo, and the re- 
sulting halo to have a much larger velocity dispersion than 
expected for an object of that mass. The massive neigh- 
bouring object will have a much higher velocity dispersion 
anyway, so will be unaffected by such effects. 

To illustrate the effect that these problems have on the 
physical properties we calculate, we show the spin distribu- 
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log A 



Figure 1. Histogram of the spin parameters from the basic FO- 
Fall halo catalogue of the Millennium Run, showing a long tail to 
high spins. The tail continues up to A ss 680, and there are over 
900 000 objects with A 0.3 (marked on the graph). The median 
spin of the distribution, A mo( j, is displayed with the uncertainty 
given by Eqn. 13. This demonstrates the need for more careful 
definition and selection of haloes. 

tion of the FOFall catalogue in Fig. 1. It shows a long tail at 
higher spins; there are 900 748 objects (6% of the catalogue) 
with spin A ^ 0.3. Fig. 2 shows the FOFall spin distribution 
as a function of halo mass, rescaled to show the fractional 
number of objects at each mass so the trend of the mass func- 
tion is removed. It shows that the high-spin tail comes from 
objects over a large range of masses, and is therefore not due 
to under-resolving groups. These presumed-anomalous high 
spins come from objects with high velocity dispersions for 
their masses, caused by situations such as those described 
above. This increases the haloes' kinetic energies, T, leading 
to large values of A. 

These features are not usually seen in other published 
spin distributions (e.g. Gardner 2001, Vitvitska et al. 2002, 
Avila-Reese et al. 2005, Tonini et al. 2006, among many oth- 
ers), partly because more advanced groupfinders and halo 
selection criteria are often used (as we do below). However, 
the fact that we can see these artifacts so clearly is because 
the Millennium Run gives us a vast number of objects, over 
a wide range of masses. 

For convenience, we shall refer to haloes suffering from 
the problems described in this subsection as "mis-defined" 
haloes, as their anomalously-high spins originate in how the 
haloes are defined by the groupfinder algorithm in relation 
to their environment, rather than any physical or numerical 
effect. 

3.2.2 A better groupfinder: The TREE haloes 

As a third definition of halo, we use the 'merger-tree' haloes 
described by Harker et al. (2006). These are the 2 = ob- 
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Figure 2. Halo spin as a function of halo mass for the FOFall 
halo catalogue. The contours indicate the relative number den- 
sity of haloes with that value of \(M^); that is, the haloes were 
binned onto a 50 X 50 grid between the maximum and minimum 
values in log Afj, and log A, and the number of haloes in each grid 
cell was normalised by the number of haloes in that mass-bin, 
thus removing the effects of the halo mass function from the plot. 
The contours are spaced logarithmically, with one contour for ev- 
ery factor of 10 0,5 in halo number density. The innermost bold 
contour (red) represents 10 _1 (i.e. a tenth of the haloes in each 
mass bin), and the outer bold contour (blue) represents 10 . 
Notice the high-spin bulge, which extends over a large range of 
halo masses. The results for SO haloes show a very similar distri- 
bution. 



jects in a catalogue of halo merger trees constructed in the 
Millennium Run (Helly et al., in preparation). These merger 
trees are similar to, but distinct from those of Springel et al. 
(2005) and Gao et al. (2005), who used different criteria for 
identifying and tracking the haloes over time. The merger- 
tree-halo catalogue used here was designed with the needs of 
the iV-body GALFORM semi-analytic models in mind (e.g. 
Helly et al. 2003, Bower et al. 2006). 

The haloes in the TREE catalogue are formed from the 
Subfind subhalo catalogue (based on the FOF haloes) in 
the following way. Each halo is initially taken to consist 
solely of its constituent subhaloes (i.e., is equivalent to the 
corresponding FOF halo, but without the "fuzz"). Next, a 
"splitting" algorithm is applied, which attempts to account 
for the linking of distinct, bound objects that often occurs 
with FOF. A subhalo can be split off from its parent halo 
if it satisfies at least one of the following conditions: 1) The 
distance between the subhalo's centre and the parent's cen- 
tre is more than twice the parent's half-mass radius; or 2) 
the subhalo has retained at least 75% of the mass it had at 
the last snapshot in which it was an independent halo. This 
second condition, which uses the merger-tree data, is based 
on the idea that if the subhalo was genuinely merging with 
the halo, it would be rapidly stripped of its mass. If however 
it was just a rapid encounter causing an artificial link, then 



its mass would be retained. If the subhalo is split off, any 
other subhaloes that reside within twice its half-mass radius 
are also split off to become part of the new halo. We shall 
refer to the catalogue of all haloes defined in this method as 
the TREEall catalogue. 

This definition of the TREE-halo catalogue alleviates 
the groupfinder-based problems somewhat. The unbound 
particles excluded from the TREE haloes have, by definition, 
higher velocities than the bound structures in the haloes. 
Their removal therefore reveals the more relaxed underlying 
haloes, with lower kinetic energies and hence A. "Interloper 
particle" contamination, in particular, is reduced by this fea- 
ture. The halo splitting algorithm also helps in many cases 
by separating objects that have been spuriously linked by 
bridges. 

The spin-mass distribution of the TREEall catalogue 
can be seen in Fig. 3; although some of the high-spin objects 
remain, the majority are now gone (there are sufficiently few 
high-spin objects remaining that they are invisible in a plot 
of P(log A) such as Fig 2). Investigating the remaining high- 
spin objects reveals them to be victims of velocity contam- 
ination from very massive neighbours. They are somewhat 
special cases however: for them not to have been rejected as 
fuzz, the contaminants must be self-bound bodies, with ^ 20 
particles each. They must also be built on density maxima 
independent of that of the host halo, in order for Subfind 
to have identified them as separate subhaloes. Furthermore, 
the interloping subhalo must be within twice the half-mass 
radius of the halo, so that the tree algorithm does not split 
it off. A consequence of this is that the resulting halo must 
consist of at least 40 particles (20 for the interloper subhalo 
and 20 for the main halo) , and this can be seen in the offset 
low-mass cutoff for high-A haloes in Fig. 3. 

3.2.3 A better halo catalogue: The QE criterion 

A relatively simple way of culling the remaining anomalous 
spin objects is to remove those that are clearly out of equi- 
librium at the moment of the snapshot. This is not quite 
the same as selecting only objects that are within a certain 
degree of true virialisation, since we don't have the neces- 
sary time-resolution to determine if the system properties 
are genuinely stationary: just as an object can appear to be 
linked to another by a bridge that may exist only fleetingly, 
so a halo could instantaneously have very similar energies to 
those of a stationary system. Therefore we will apply a cut 
in the instantaneous 'virial ratio' of halo energies, 2T/U + 1, 
and describe the objects that meet this criterion as haloes in 
a quasi- equilibrium (QE) state. This name avoids implying 
the zero time-derivative necessary for the true virial ratio. 

The question of where to make the cut in 'virialisation' 
(i.e. applying a QE limit) is a difficult one, because the deci- 
sion will always be somewhat arbitrary. Since it is desirable 
to minimise the effect of such arbitrariness, we concentrate 
on the effect of applying a QE limit to the TREE haloes, 
since the merger-tree algorithm has already removed many 
of the mis-defined haloes. 

The value of 2T/U + 1 for the TREE haloes is plotted 
against their mass in Fig. 4. We applied a QE cut of the 
form 

2T 

-Q<— + 1^Q (12) 
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Figure 3. Spin as a function of halo mass for haloes in the 
TREEall catalogue. The contouring is as in Fig. 2, i.e. in equal 
logarithmic steps of 10 0,5 , normalised to remove the mass func- 
tion. In this plot, the inner bold contour (red) represents 10" 1 
of the haloes at each mass, and the outer bold contour (blue) 
represents 10~ 5 of the haloes at each mass. The merger-tree halo 
definition has moved many of the high-spin haloes visible in Fig. 
2 down into the main body of the distribution. 



to the TREEall catalogue, examining the effect of a wide 
range of Q-values on the P(logA) and A(Mh) distributions. 
Because of the relatively small numbers of objects with 
anomalously high spins, we find that the QE cut makes 
negligible difference to the shape of the spin distribution, 
P(logA). A very small value, Q < 0.3, will act to shift the 
median spin lower by a few percent, due to the mass de- 
pendence seen in Fig 4. Through a detailed examination of 
the A(Mh) distribution, we find that a value of Q between 
0.4 and 0.6 gives a good balance between removing the mis- 
defined haloes and reducing the overall sample size (adding 
noise and biasing it by mass) . Higher values of Q allow some 
mis-defined haloes to creep in, with a significant impact for 
Q > 1.0. We will use a value of Q = 0.5 for our cleaned 
halo catalogues. This cut is shown in the horizontal dashed 
lines of Fig. 4, and the resulting A(Afh) distribution is shown 
in Fig. 5. Applying this form of virialisation cut on the halo 
catalogue provides a useful tool with which to remove haloes 
with anomalous spins caused by mis-defined haloes. 



3.2.4 A better halo catalogue: Numerical effects 

A second peculiarity of the spin distributions is visible 
in Figs. 2, 3 and 5 (for the FOFall, TREEall and quasi- 
equilibrium TREE halo catalogues respectively): an upturn 
in the spin distribution at low masses. This can be seen 
clearer in the variation of the median spin over mass bins 
A m ed(Mh), plotted for the FOF haloes in Fig. 6. This effect 
is unrelated to the velocity contamination problems of the 
mis-defined haloes, and instead comes from the mass resolu- 
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Figure 4. Contour plot of the instantaneous 'virial ratio', 2T/U + 
1, against halo mass for TREE haloes. A virialised object has a 
value around zero, and a gravitationally bound object has value 
> — 1. The tail at low values (large T) extends down to 2T/U+1 ~ 
-960; there are 3733 objects with 2T/U + 1 < -1. The dashed 
lines show the QE limit of Q = 0.5, and the lower particle-number 
limit of N p = 300. The contouring is as in Fig. 2, i.e. relative halo 
number density in equal logarithmic steps of size 10 ' 5 . The inner 
bold contour (red) represents 10 _1 of the haloes at each mass, 
and the outer bold contour (blue) represents 10 — 4 of the haloes 
at each mass. The plots for FOF and SO groups are similar to 
this. 



tion of the simulation affecting the angular momenta. This 
effect has been seen before, for example by Reed et al. (2005) 
in the context of subhaloes. To understand the cause of this 
effect, consider a continuous object with angular momentum 
Jtrue- If we construct a realisation of this object using a sam- 
ple of N discrete particles, the resulting angular momentum 
can be modelled as the vector sum of the 'true' angular mo- 
mentum (from the continuous object) with a noise vector 
oriented in a random direction: J — Jtrue + Jnoise- This 
will act to push the measured magnitude J up above Jtrue 
because the random direction of Jnoise will mean it reaches 
outside the sphere of radius j Jtrue | more than 50% of the 
time. Therefore, the random noise inherent in using discrete 
particles to sample a near-continuous object such as a dark 
matter halo would act to bias J upwards, and A along with 
it. 

Modern A-body codes such as the L-GADGET-II code 
used for the Millennium Run are very good at conserving 
quantities such as energy and angular momentum, so for 
a well-resolved object there is negligible inaccuracy aris- 
ing from particle discreteness. For a less well-resolved ob- 
ject however, the effect can nevertheless be relevant, even 
though the angular momentum of the particles making up 
the halo has been well conserved. Discreteness mainly af- 
fects the outer parts of a halo, making the effective surface 
more jagged than that of the continuous object it represents. 
We also expect the outer parts to harbour most of the an- 
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Figure 5. Spin parameter against halo mass for TREE haloes 
with a quasi-equilibrium (QE) limit of 0.5 applied; this can be 
compared with Figs 2 and 3. The same contouring is used, i.e. 
relative halo number density in equal logarithmic steps of size 
10 ' 5 . The inner bold contour (red) represents 10 — 1 of the haloes 
at each mass, and the outer bold contour (blue) represents 10~ 4 
of the haloes at each mass. The QE-limit has removed the vast 
majority of the high-spin haloes seen in Fig. 3. 




log N p 

Figure 6. Median spins for all FOF haloes for the Millennium 
Run (thick red line) and the milli_lowres run (thin blue line). 
The most massive bin contains 2000 haloes, and the remaining 
bins are logarithmically spaced, with one bin every factor of 2 
in halo mass (except for the leftmost bin which extends down 
to the cutoff limit of 20 particles.) Horizontal bars mark the bin 
widths, and vertical error bars give a measure of the precision of 
the median according to Eqn. 13. 



gular momentum. For lower-mass haloes, a greater fraction 
of their particles make up these 'surface' regions, so this 
problem has a greater effect; the inclusion of a single par- 
ticle can make a significant contribution to A. Hence, the 
spin parameter rises for haloes with fewer particles because 
the discreteness of the haloes' surface layers adds noise to 
their 'true' angular momenta. (This is not the same effect 
as discussed in Shaw et al. (2006); there, a surface pressure 
term is added to the virial theorem to account for their halo 
truncation at R v i r .) 

The importance of the noise contribution to A can 
be examined by determining the spin distribution of the 
same simulation (same code and same corresponding ini- 
tial condition waves) but ran at a different resolution. 
We performed a lower resolution resimulation of the Mil- 
lennium Run, with 2160 3 /64 = 540 3 particles (so their 



mass is m p .i ow = 64m 



p,Millen ) , 



which we will refer to as 



milli_lowres. The FOF and Subfind algorithms were im- 
plemented on milli_lowres in the same way as in the Mil- 
lennium Run itself. (Although we do not have merger-tree 
data for milli_lowres, and hence cannot construct a TREE 
halo catalogue, the Millennium Run results show that the 
same effect is seen in FOF, SO and TREE haloes.) Fig. 
6 shows the median spin A me d(A r p ) for FOF haloes in the 
Millennium Run and milli_lowres, with vertical error bars 
showing a Gaussian-like estimate of the precision of the me- 
dian given by: 

(A84 — Al6j 



fmcd 



where Ai is the ith percentile of A 



(13) 



i - 16% = 68%, the 



amount of data within ±lcr of a Gaussian peak) and iVhaio is 
the number of haloes in that mass bin. Note that the spread 
in A is much greater than the uncertainty in A me d; compare 
the error bars in Fig. 6 with the data shown in Fig. 2. 

Although the two curves show the same qualitative be- 
haviour (the low-iVp upturn), there is a vertical shift be- 
tween them. This is due to milli_lowres containing fewer 
of the mis-defined objects described in Section 3.2.1 than 
the Millennium Run itself. Decreasing the resolution effec- 
tively smooths the density field, so that small objects with 
more massive neighbours can disappear completely, whereas 
a more isolated object of the same mass may still survive 
(although containing fewer particles). This means that the 
"real" objects are retained (and there are still many under- 
resolved ones causing the upturn in A me d(A'p)), but there 
is a reduction in the number of mis-defined objects. The 
milli_lowres results seem to confirm the dominance of nu- 
merical effects at low-iVp, above which the effect of noise is 
negligible. 

If we apply the quasi-equilibrium cut described in sec- 
tion 3.2.3 to remove the mis-defined objects, we can examine 
the effect of discreteness on the median spins of just the 'real' 
haloes. Fig. 7 shows A mo( j(A'p) for QE-selected FOF haloes, 
compared with QE-selected FOF haloes from milli_lowres. 
In contrast to Fig. 6, the two lines now lie on top of each 
other, exhibiting the same upturn in spin for haloes with 
the same number of particles. This demonstrates that the 
upturn is purely a numerical effect. We can exclude haloes 
that appear to be dominated by this effect by fixing a limit 
of A r p > 300 on the halo catalogue. 
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Figure 7. Median spins against number of particles in haloes, for 
FOF haloes with the QE cut applied. The thin blue line shows the 
milli_lowres haloes and the thick red line shows those from the 
Millennium Run itself. The error bars show the uncertainty on 
the median, using Eqn. 13. The two lines show an identical trend 
at low ./Vp, demonstrating that the upturn in A mo( j is indeed a 
numerical effect, affecting the spins of haloes containing fewer 
than about 300 particles. 



Fig. 4 shows how the QE and JV P cuts we use relate 
to one another for the TREE haloes. The low-iV p cut has 
a significantly stronger effect on the halo catalogue. Ap- 
plying the QE cut on its own reduces the population of 
the TREEall catalogue by only 0.3%; applying the \aw-N p 
limit as well means removing in total 91.17% of the origi- 
nal TREEall haloes (see Table 2) . We refer to the resulting 
cleaned TREE-halo catalogue as TREEclean, and it is these 
haloes whose properties we shall be examining in detail. 
In some cases, for completeness, we shall compare the re- 
sults with those from the FOFclean and SOclean catalogues. 
These are cleaned using the same Q = 0.5 and N p ^ 300 
cuts as the TREEclean catalogue. 

It is important to note that the criteria we have adopted 
- the QE cut and low-jV p limit - are those appropriate to the 
quantities of interest in this work. For example, Fausti-Neto 
et al (in prep.) use the substructure parameter S = Ar/i? v j r , 
where Ar is the distance between the potential minimum of 
the halo and its centre of mass. Their final criteria for se- 
lecting haloes in the Millennium simulation are S < 0.07 
and + 1 > -0.25. They are concerned with fitting den- 
sity profiles to haloes, and using the substructure parameter 
allows them to remove haloes with a large fraction of mass 
in substructures which would otherwise contaminate their 
results. Similarly, Maccio et al. (2006) define an "offset pa- 
rameter" as x H = Ar/i? v i r , where Ar is measured from 
the most-bound particle rather than the potential minimum. 
They use this alongside the rms of their density profile fits, 
to assess the quality of their halo catalogues for estimates 
of halo concentration. Although we have examined the sub- 



structure/offset parameter and how it affects A, we find that 
it is not useful in removing the mis-defined haloes, or those 
whose spins are dominated by the numerical effects discussed 
above. 

Having successfully implemented an appropriate 
groupfinder and cleaned the resulting halo catalogues, we 
can now proceed to examine their spin properties. The FOF 
and TREE halo catalogues, including some halo properties 
and semi-analytic galaxy properties, are publicly available 
online 1 (Lemson et al. 2006). 



4 RESULTS 

4.1 The form of the spin distribution 

The median spin of the TREEclean halo catalogue is A me d = 
0.0381. The distribution of halo spins about the median, 
-P(A), has been often fitted with a lognormal function (e.g. 
van den Bosch 1998; Gardner 2001; Bailin & Steinmetz 
2005), i.e. a Gaussian in log A: 



P(logA) = 



(Tig 



exp 



log 2 (A/Ao) 



2< 



(14) 



While this fitting function has proved adequate for small 
numbers of objects, we find that for the > 10 6 haloes 
in the Millennium simulation, deviations from a Gaussian 
are clear and significant. The spin distribution drops faster 
than a Gaussian at high spins, and slower than a Gaus- 
sian at low spins. The best fit to the TREEclean catalogue 
is shown in Fig. 8, which fits Eqn. 14 with peak location 
A = 0.03687 ±0.000016 and width er lg = 0.2216 ± 0.00012. 2 
The corresponding lognormal function of A has the same 
peak, and a width of a — ln(10)o"i g . The fit has a reduced- 
X 2 of 40.46. 

Part of the reason why a lognormal is such a poor fit 
is that this function strongly avoids very low spin values, 
whereas the real distribution, based as it is on the three- 
dimensional vector j , does not. The longer tail at low- A is 
primarily due to the distribution of j being smooth and 
isotropic about j = 0, implying 3 that P(logA) oc A 3 . 

We have found that the following function provides a 
better description of the data: 



(15) 



For the normalised spin distribution, we can express A in 
terms of the other free parameters, a and Ao (the peak lo- 
cation): 







' A \ 3/a " 


(£) exp 


-a 





A = 3 In 10 



r(a) 



(16) 



1 http:/ /www. mpa-garching.mpg.de/millennium/ 

2 Throughout this paper, the quoted uncertainties on best-fitting 
parameters are given by the square root of the diagonal of the 
covariance matrix for that fit. 

3 Maccio et al. (2006) claimed the low-A tail is due to the higher 
uncertainty in A at low values. However, by varying the minimum 
N p for the halo catalogue, and hence the uncertainty in A, we 
found that the low-A side of the distribution consistently drops 
off slower than the high-A end, confirming that this shape is not 
primarily due to uncertainties. 
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where the gamma function T(pt) = (a — 1)!. The best fit to 
the data is shown in Fig. 9, and has parameters: 



A„ = 0.04326 ± 0.000020 



a = 2.509 ± 0.0033 



with a much-improved reduced- x 2 of 2.58. 

We have examined whether the deviations from lognor- 
mal depends on our choice of the quasi-equilibrium param- 
eter Q when cleaning the halo catalogue. We found that 
the fit remains good over a wide range of Q. We have also 
found the best fit of Eqn. 15 to the FOFclean and SOclean 
catalogues. The results for SOclean haloes are remarkably 
similar to those for the TREEclean haloes, with a reduced- 
X 2 of 3.10: 



Ao = 0.04174 ± 0.000022 



a = 2.540 ± 0.0036 



The median spin of the distribution is A me d = 0.0367. The 
haloes in the FOFclean catalogue, which have a median spin 
of A me d = 0.04288, are not as well fitted by Eqn. 15. The 
reduced-;^ 2 is 15.0 and the parameter values are: 



Ao = 0.04929 ± 0.000027 



a = 3.220 ± 0.0046 



This is, in fact, slightly worse than the best-fitting lognormal 
(Eqn. 14), which yields a reduced-x 2 of 12.5, with a peak 
location A = 0.04222 ± 0.000022 and width <7i g = 0.2611 ± 
0.00016. 

These tests show that the distribution of A depends on 
the careful definition and selection of dark matter haloes. 
However, the fact that the distribution is non-lognormal is 
not a consequence of the particular choice of groupfinder or 
selection criteria used here - the form of the distribution is 
not peculiar to the TREEclean catalogue. 
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Figure 8. The normalised spin distribution (red histogram) for 
merger-tree haloes with N p > 300 and filtered using the QE crite- 
rion (the TREEclean catalogue). The best-fitting single Gaussian 
function is plotted as the smooth blue curve. The Gaussian drops 
too quickly at low spins, and too slowly at higher spins. To try 
to minimise these effects, the best fit has a peak location that is 
shifted away from that of the histogram data. The median spin 
of the distribution is also displayed. 



4.2 Correlation of spin and halo mass 

The variation of median spin parameter with halo mass, for 
the cleaned catalogues from the three groupfinders is shown 
in Fig. 10. ft is interesting to note that the FOF haloes 
exhibit an upturn in spin for objects more massive than the 
low-iVp cut, an effect that is not present in the TREE or 
SO haloes. This can be attributed to the outer parts of the 
FOF haloes consisting mainly of unbound 'fuzz' particles. 
These will usually have higher velocities, which act to innate 
the spin. These particles are not part of the TREE haloes, 
and most will be shaved off in SO haloes too. The SO and 
TREE haloes show a shallow downwards trend of A mc d up 
to M h ~ 10 13 /i _1 Mq ~ 12000m p and a rapid decline at 
larger masses. 

We fit a cubic polynomial to the TREEclean median 
spin data, 

log A mcd = ax 3 + fix 2 + 72; + S (17) 

where x — log Afh/( /i _1 Mq). The best-fitting values of 
these parameters are: 

a = (-8.6 ± 1.4) x 10~ 3 

fi = (3.2 ± 0.54) x lO -1 

7 = -4.1 ±0.68 

5 = 15.7 ±2.8 

with a reduced-^ 2 of 0.44. 

While the trend of A mc( j with mass is real, it is impor- 
tant to note that it is a small effect; the scatter around this 




log A 



Figure 9. Histogram (red) of the normalised log A distribution of 
TREEclean haloes, as in Fig. 8. The smooth solid curve (blue) is 
the best fit to the data using Eqn 15. The peak location of the fit 
(Ao) and the median of the data (A med ) are marked with arrows. 
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Figure 10. The median spins of halo mass bins, for TREEclcan 
haloes (bold blue points), SOclean haloes (medium green line) 
and FOFclean haloes (thin red line). Data have also been calcu- 
lated below the 300-particlc low-mass limit (marked with a dashed 
line). The thick black curve is the best-fitting cubic polynomial to 
the TREEclean data. The vertical error bars are given by Eqn 13. 
The mass-binning scheme is similar to Fig. 6, but with bins every 
factor of \/2 in mass. Only the widths of the most massive two 
bins are marked, for clarity. 



Figure 11. Axis ratios p and q for the cleaned TREE, FOF and 
SO halo catalogues. Prolate objects have p = 1, oblate objects 
have 5 = 1, and spherical objects have p = q = 1. The SO haloes 
are more spherical than the other two types. TREE and FOF 
objects exhibit a range of shapes, and all three catalogues show 
a preference for prolatcncss over oblateness. The contouring is by 
halo numbers, in equal logarithmic steps of 10 5 . In each plot, 
the red inner bold contour represents 10 3 haloes, and the blue 
outer bold contour represents 10 1,5 haloes. 



median is large (compare with Fig. 5, which shows data from 
the same haloes but on a log A scale). This is in qualitative 
agreement with previous results (e.g. Cole & Lacey 1996), 
but because of its weak nature, this trend has often not been 
visible (e.g. Warren et al. 1992; Lemson & Kauffmann 1999; 
Maccio et al. 2006). 

4.3 The halo shape distribution 

The shapes of the haloes are described by the axes, a ^ 
b ^ c, of the ellipsoid derived from the inertia tensor, as 
described in Section 2.2. Fig. 11 shows the distributions 
of p = c/b and q = b/a for the three cleaned halo cata- 
logues. The haloes are generally triaxial, but they have a 
range of shapes, with a slight preference for prolateness over 
oblateness. The distribution agrees qualitatively with previ- 
ous work such as that by Frenk et al. (1988), Warren et al. 
(1992), Cole & Lacey (1996), Faltenbacher et al. (2002) and 
Bailin & Steinmetz (2005). Unsurprisingly, SO haloes are 
more spherical than FOF or TREE haloes. FOF haloes show 
a much broader distribution of shapes (and a stronger pref- 
erence for prolateness) than SO or TREE haloes. 

Fig. 12 shows how the median shape of haloes in the 
TREEclean catalogue changes with halo mass, using the 
minor-to-major axis ratio s and the triaxiality parameter 
T = (a 2 — b 2 )/(a 2 — c 2 ). More massive haloes tend be 
less spherical and more prolate. Again, this is in qualita- 
tive agreement with previous results, such as those of War- 
ren et al. (1992), Bullock (2002), Kasun & Evrard (2005), 



Shaw et al. (2006), Gottlober & Turchaninov (2006), All- 
good et al. (2006), and Maccio et al. (2006). This is also 
what one might expect in a hierarchical formation model in 
which haloes tend to form by matter collapsing along fila- 
ments, leading to prolateness, rather than onto sheets which 
would lead to oblateness. Furthermore, the more massive 
haloes form later, and have had less time to relax into more 
spherical configurations. Since we have deliberately tried to 
select the more relaxed objects, the remaining trend we see 
here is weak. Furthermore, although the medians follow well- 
defined trends, the spread of the distribution in halo shapes 
covers virtually the entire range in both s and T, as can be 
seen from the percentile bars on the graph. 

The two graphs in Fig. 12 both show a change in be- 
haviour not seen in previous work, around the N p = 300 
limit. Resolution tests similar to those described in section 
3.2.4 were carried out using the milli_lowres simulation, to 
assess whether this change in behaviour was indeed a numer- 
ical effect similar to that seen in halo spins (in Figs. 6, 7 and 
10). The results showed that these halo shape parameters do 
also require N p > 300, reinforcing our previous choice. In- 
deed, one would expect haloes whose spins are affected by 
particle discreteness (i.e. with N p < 300) to be less spherical 
and more stringy (prolate). 

We fit a broken line to both shape parameters for the 
TREEclean catalogue, of the form: 

J mix + ci x < x . 

Vmedlx) = < (18) 
W ^ m,2X + C2 X > x 
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Figure 12. The median, in bins of halo mass, of the axis ratio s 
(top) and the triaxiality parameter T (bottom), for TREEclcan 
haloes, using the same mass binning scheme as in Fig. 6. The 
low-N p limit of 300 particles is marked with a dashed line, and 
data points have also been plotted below this limit). The medians 
for FOFclean and SOclcan haloes follow the same behaviour as 
the TREEclean haloes, but with SO objects being more spheri- 
cal and FOF objects less; FOF haloes are more prolate and SO 
haloes show a weaker preference for prolateness. Error bars on 
the medians (following Eqn. 13) are plotted, but most are van- 
ishingly small. The whisker bars are percentiles, at the equivalent 
of lcr (68% of haloes, boxes) and 2<r (95% of haloes, bars); from 
the bottom to the top of each graph, these show where approx- 
imately 2.5%, 16%, 84%, and 97.5% of the haloes have s (or T) 
below these values. The thick blue lines show broken-line fits to 
the data; see the text for details. The dotted red line is the fit of 
Allgood et al. 2006, who used a different definition of a halo (see 
text) . 



where x = logMh/( /i -1 Mq). We fit with mi, m?, xo and C2 
as free parameters, with Ci = c 2 + (m 2 — mi)xo. The fitted 
parameters for s me d are: 



mi,, = (9.2 ±0.87) x 10 , Ci,. 
m 2 ,s = (-6.6 ±0.12) x 10" 2 , c 2 „ 
x 0s = 12.27 ±0.012 



0.56 ±0.015, 
1.48 ±0.015, 



with a reduced- x of 29.9. The fitted parameters for 7Jn et i 
are: 



mi,T = (-1.6 ±0.18) x 10" 2 , 
m 2 , T = (7.2 ±0.24) x 10~ 2 , 
x ,r = 12.28 ±0.021 



ci, r = 0.82 ±0.031, 
c 2 ,T = -0.25 ±0.029, 
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Figure 13. The main plot shows spin versus axis ratio s = c/a 
for TREEclean haloes; contouring is similar to Fig. 2, showing the 
number density of haloes, normalised by the number of haloes in 
each s-bin. The red inner bold contour represents 10 _1 of the 
haloes in each s-bin, the blue outer bold contour represents 10 — 3 
of the haloes in each s-bin, and the contours are spaced in equal 
logarithmic steps of 10 0,5 . The upper plot is a histogram of s for 
TREEclean haloes, effectively showing the function by which the 
contour plot has been normalised. 



with a reduced-^ 2 of 4.27. The two mass breakpoints xq jS 
and xo t r agree within their uncertainties. 

Allgood et al. (2006) fit a power-law to s me d(Mh). This 
is plotted in Fig. 12, and indicates that that their haloes are 
significantly less spherical than ours. This is largely a result 
of different group definitions; although not plotted, we find 
that our SOclean and FOFclean catalogues differ from the 
TREEclean results by a similar amount. A power-law of the 
type used by Allgood et al. (2006) would not be a good fit 
to the data presented here which have a definite change in 
slope towards lower halo masses. 

The overall distributions of s and T are shown in the 
upper plots of Figs. 13 and 15. These agree well with distri- 
butions seen previous work, e.g. Bailin & Steinmetz (2005) 
and Shaw et al. (2006). The drop-off in halo sphericity below 
about s ~ 0.3 can be explained by considering how flatter 
haloes would puff up due to bending instabilities (Merritt & 
Sellwood 1994). 



4.4 Spin and shape parameters 

The relationship between spin parameter and halo shape is 
illustrated in Figs 13, 14 and 15. Fig. 14 emphasises the 
trend visible in Fig. 13 by plotting the median spin parame- 
ter as a function of s for different mass bins. There is a clear 
trend for more spherical haloes to exhibit less coherent ro- 
tation. Although this trend is in the sense one might naively 
expect, the haloes, in fact, do not have very high spin, and 
are not rotationally supported. The origin of this trend is 
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Figure 14. The median spin parameter in bins of axis ratio 
s = c/a, for haloes from the TREEclean catalogue. The heavy 
black line is for all haloes with more than 300 particles (the full 
TREEclean sample); the other lines show the trends for different 
halo mass bins. The horizontal bars show the widths of the first 
and last bins for each line. The trend is for more spherical haloes 
to have less coherent rotation in the median. This trend becomes 
very steep for the most aspherical haloes (although these are not 
present at higher masses because of the rapid drop in the halo 
mass function). 

likely to lie instead in the effects of the tidal torques experi- 
enced by the haloes during their early phases of formation. 

As seen previously in Fig. 10, the least massive objects 
have the most extreme spins, in the median. Fig. 14 shows 
that the higher spin objects are also less spherical. The 
haloes which are closest to spherical have a spin parame- 
ter that is independent of halo mass, and has A mcd ~ 0.033. 
(This does not apply to the most massive haloes, however, 
since their population lacks the more spherical objects) Fur- 
thermore, the median spins for the more massive haloes are 
independent of shape and have A mc( j « 0.032. 

In contrast to the variation of A mc( j with s, Fig. 15 shows 
that there is only a very weak trend of spin with halo triaxi- 
ality. Over the entire range of triaxiality, each T-bin contains 
a very similar fraction of haloes at each value of log A. 

4.5 Spin-shape alignment 

Fig. 16 shows the angle between the angular momentum 
vector and the three shape axis vectors, e.g.: 

cos6> a = |j'-a/j| (19) 

for alignment with the semimajor axis given by the unit vec- 
tor a. Note that this definition does not distinguish between 
j lying parallel or antiparallel to the axis vectors. 

Most haloes have their spin axis well aligned with their 
minor axis, and lying perpendicular to their major axis. 
However, the distribution of alignments with respect to all 
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Figure 15. The main plot shows spin versus triaxiality parameter 
T for TREEclean haloes; contouring is as for Fig. 13, i.e. in equal 
logarithmic steps of 10 ' 5 haloes per T-bin. The red inner bold 
contour represents 10 -1 of the haloes in each T-bin, and the 
blue outer bold contour represents 10 -3 of the haloes in each T- 
bin. The upper plot is a histogram of T for TREEclean haloes, 
effectively showing the function by which the contour plot has 
been normalised. 
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Figure 16. Normalised histograms of the cosines of the angle 
between the angular momentum vector and the major (thick red 
line), intermediate (medium green line) and minor (thin blue line) 
axes of haloes in the TREEclean catalogue, as defined in the text. 
A random distribution would be a flat line at P(cos 9) = 1. 
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Figure 17. Normalised histograms of the cosine of the angle 
between the specific angular momentum vector and the major 
(thick red line), intermediate (medium green line) and minor (thin 
blue line) axes of the haloes. The TREEclean catalogue has been 
cut at the values of the triaxiality parameter shown. The most 
spherical haloes (s > 0.8) have been removed from each sample. 
The 'very prolate' selection only contains 1360 groups, giving rise 
to a noisier histogram. 



three axes is fairly broad. This agrees with previous results 
(e.g. Warren et al. 1992, Bailin & Steinmetz 2005, Allgood 
et al. 2006 and Shaw et al. 2006). 

Extremely oblate objects have a degeneracy between 
the major and intermediate axes ( £ ~ 1), so there is an equal 
probability for the angular momentum vector to subtend a 
given angle with either axis This can be seen in the top-left 
panel of Fig. 17, which shows the alignment distribution for 
the more oblate haloes, i.e. those with T < 0.2 (and s < 0.8 
in order to avoid an a ~ b ~ c degeneracy). Since most 
haloes have j aligned with c, j has a preference for being at 
right angles to the two large axes. 

For extremely prolate haloes, the degeneracy is between 
minor and intermediate axes (| ~ 1). In this case, the dis- 
tribution of the alignment of the angular momentum vector 
with either of these axes is similar only for the tail of ex- 
tremely prolate haloes, T > 0.99 (bottom-right panel). For 
T > 0.8, only a small number of haloes have their spin axis 
aligned with the intermediate axis (top-right panel); like the 
bulk of the population as a whole, most of these haloes ro- 
tate around their smallest axis. 

The reason for the distinction between the major and 
intermediate axes present in Fig. 16 can now be seen to be 
a combination of various effects. The preference for prolate- 
ness over oblateness means that there is a tendency for the 
intermediate axis to be more similar to the minor than to 
the major axis. This increases the probability of j being 
aligned to the b axis. However, since j can be aligned with 
only one axis (and when this happens it must be perpen- 
dicular to the other two), the preference for the spin axis to 



be the minor axis outweighs the preference for prolateness, 
and this results in the slight excess probability for j to be 
perpendicular to the intermediate axis seen in the figure. 



4.6 Spin, shape and Halo clustering 

In this section, we investigate whether halo spin and shape 
has an environmental dependence. We quantify environment 
by means of the haloes' two-point correlation function, f (r), 
and we explore whether the clustering amplitude differs for 
haloes with different spin. The three halo catalogues were di- 
vided into four bins in mass, and the haloes in each mass bin 
were then divided about the median spin for that mass; f (r) 
was calculated for each set of haloes. The results are plotted 
in Figs 18 and 19 for the TREEclean and FOFclean cata- 
logues respectively. (The results for SOclean haloes are very 
similar to the FOF results.) The lines in Fig. 18 show power- 
law fits over the limited range of pair separations shown, 

The fitted values of ro and 7 are given in Table 3. 

For the higher mass bins, the results from the differ- 
ent types of haloes are similar: higher spin haloes are more 
strongly clustered than lower spin haloes. This could be be- 
cause objects evolving in denser, more clustered environ- 
ments are more likely to experience stronger tidal forces, 
leading to more coherent rotation. These results are consis- 
tent with the tentative earlier work of Barnes & Efstathiou 
(1987), as well as the mark correlation function analysis of 
Faltenbacher et al. (2002) 

The difference in clustering strength between high and 
low spin haloes decreases with halo mass. For the least mas- 
sive TREE haloes, M = 10 n ' 4 -10 12 h -1 MQ, there is virtu- 
ally no difference in the correlation functions of the fast and 
slow rotators. 

Note that haloes from the 'cleaner' groupfinder (TREE) 
exhibit stronger clustering at all scales, for all but the high- 
est masses. This is due to the fact that the splitting algo- 
rithm applied during the construction of the merger trees 
results in a greater number of close halo neighbours (albeit 
with reduced masses) compared to the corresponding FOF 
haloes. 

We have performed a similar analysis to that presented 
in Figs. 18 and 19 for the halo sphericity parameter, s. We 
found analogous results, in that the more spherical haloes 
are more clustered than the less spherical haloes. 

To examine the effect of halo spin and shape on cluster- 
ing in more detail, we consider the bias parameter, 6, which 
describes how much more or less clustered a set of haloes 
is relative to the underlying dark matter distribution. We 
examine how the bias varies for haloes with different values 
of A or s, at a fixed range of mass. Similar analyses have re- 
cently been performed by Gao et al. (2005), Wechsler et al. 
(2005) and Wetzel et al. (2006), who examined the effect 
of halo formation time, concentration, substructure content 
and time since last major merger on the bias. 

The bias parameter is related to the correlation function 
through: 

£hh(r|A/ h , A) = b 2 (r\M b , A)£ mm (r) (21) 
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Mass bin 


A < A 


tncd 


A > A 


med 


(fc _1 Mg) 


7 




7 




10 11.4_ 10 12 


-1.553 ±0.0013 


3.670 ± 0.0070 


-1.489 ± 0.0012 


3.671 ± 0.0073 


10 12 -10 13 


-1.591 ± 0.0031 


4.22 ± 0.020 


-1.537 ± 0.0026 


4.82 ± 0.021 


10 13 -10 14 


-1.64 ± 0.021 


6.6 ± 0.24 


-1.71 ± 0.016 


8.3 ±0.23 


10 14_ 10 16 


-1.6 ±0.19 


13 ± 5.8 


-1.7 ± 0.12 


19 ±4.8 



Table 3. Parameters for the power-law £(r) = (r/r ) 7 from fitting to the eight two-point correlation functions for TREEclean haloes. 
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Figure 18. Two-point correlation function, £(r), for TREEclean haloes in four mass bins. The lines represent power-law fits over the 
r-bins shown. The thick lines and red points are for haloes with A ^ A mct j for that mass bin; the thin lines and blue points are for 
haloes with A < A mec j. The data are noisy in the higher mass bins which contain fewer haloes. The error bars are Poisson errors, i.e. the 
square-root of the number of pairs in each r-bin, divided by the mean number of pairs per r-bin. 



where ^^(rlM^, A) denotes the halo-halo correlation func- 
tion for haloes in a given range of mass and spin (in this 
example), and £ mm (r) is the dark matter correlation func- 
tion. We compute the bias parameter as a function of mass, 
6(Mh), using a similar method to that of Gao et al. (2005). 
Specifically, we compute £hh(?"|A^h, A) in four r-bins in the 



range 6 ^ r ^ 25/t _1 Mpc, equally spaced in logr. The bias 
parameter at each mass is then found as the normalisation 
constant that minimises: 

2 _ / £hh(rQ - ^mmtr.) \ 2 / 00 s 

Z — 1 ' 
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Figure 19. As Fig. 18, but for the FOFclean haloes. No fits were made to these data, and the lines merely join the points. For this 
catalogue, the higher spin haloes (thick red lines) are consistently and significantly more strongly clustered that the lower spin haloes 
(thin blue lines). The results for the SOclean haloes are very similar. 



where the <7hh(?"i) are the Poisson errors on £hh(Vi), and the 
sum is over all (> 1) r-bins where £hh(n) > and £hh(Vi) 
was computed using at least 100 objects. This procedure is 
performed first for all the TREEclean and FOFclean haloes. 
It is then repeated for the haloes in the upper and lower 
20th percentiles of the A and s distributions, from both cat- 
alogues. 



The results are shown in Fig. 20. This shows that 
higher spin haloes are more clustered than the average, 
and the lower spin haloes are less clustered. This trend is 
largest at higher masses, reaching a factor of « 2.2 be- 
tween the high-spin and low-spin bias at the largest mass, 
~2x W 14 h~ 1 M Q . 



5 CONCLUSIONS 

The huge size and high resolution of the Millennium simu- 
lation (Springel et al. 2005) makes it possible to determine 
the properties of dark matter haloes in the ACDM cosmol- 
ogy with unprecedented statistical power. In this paper, we 
have concentrated on the spins and shapes of dark matter 
haloes, ranging in mass from those of dwarf galaxies to those 
of rich clusters. We have investigated the distribution of the 
spin parameter, A, its dependence on mass, the distribution 
of shapes, the relationship between shape and spin, and the 
environmental dependence of spin and shape. We provide 
accurate fitting formulae for several quantities of interest. 

While many of the properties we have investigated here 
have been studied in earlier simulations going back over 
twenty years, a novel aspect of our work is the analysis 
and comparison of haloes identified in different ways. Along- 
side the traditional "friends-of-friends" (FOF) algorithm of 



Spin of dark matter haloes 17 



12 13 14 




'— 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 " 

-10 12 

log( M h /M» ) 



Figure 20. Bias parameter, fe(Mh), for the halo populations indicated in the legend. In each plot, the central dotted line marks the bias 
for the entire population at that mass. The thick red line is the bias for haloes in the upper 20 per cent of the distributions of A (upper 
two plots) and s (lower two plots), and the thin blue line is for haloes in the lower 20 per cent of the distributions. The error bars give 
the Ax 2 = 1 confidence interval. The lines stop either when the correlation function f(r) for at least 3 of the 4 radial bins is non-positive, 
or if it was made using fewer than 100 objects. 
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Davis et al. (1985) and the "spherical overdensity" (SO) 
algorithm of Lacey & Cole (1994), we have introduced a 
new halo definition, the TREE haloes, which are perhaps 
the most appropriate when carrying out comparisons of the 
simulation results with galaxy and cluster data. The TREE 
haloes are denned as branches of the halo merger trees, in 
which special care has been taken to identify physical haloes 
by separating objects that are artificially and transiently 
linked together. Each groupfinder results in more than 15 
million haloes at z = 0. 

The TREE halo catalogue was further 'cleaned' in 
two ways. Firstly, to remove any remaining spurious ob- 
jects, we applied a cut in the instantaneous virial ratio, 
(-0.5 «S ^ + 1 ^ 0.5). Secondly, to remove objects whose 
angular momentum is biased due to particle discreteness, 
we considered only haloes with more than 300 particles, as 
indicated by convergence tests. Our final cleaned halo cata- 
logues consist of > 10 6 haloes at z = 0. 

We find that the distribution of the dimensionless spin 
parameter, -P(A), is poorly fit by a lognormal when this 
many objects are considered. The function given by Eqn. 15 
provides a much better description of the data. Although the 
distribution of A(Afh) is fairly broad, there is a clear trend of 
the median spin with halo mass, with more massive haloes 
spinning more slowly. However, the strength and shape of 
the trend is significantly different for different halo defini- 
tions. The cubic polynomial of Eqn. 17 provides a very good 
fit to the median spin of the TREEclean halo catalogue, over 
a factor of ~ 10 3 in halo mass. 

We analysed the shapes of the haloes, and found, as in 
previous studies, that there is a broad distribution of shapes 
with a slight preference for prolateness over oblateness. More 
massive haloes are less spherical and more prolate in the me- 
dian, although the data span a large fraction of the available 
shape parameter space. We fit broken lines to the trends 
with log A/h of the median sphericity axis ratio s = c/a and 
the median triaxiality parameter T (Eqn. 18). Both these 
quantities exhibit a change of behaviour at a galactic mass 
scale, Aih ~ 2 x 10 12 /i _1 Mq, where the gradient of the fit 
changes sign, with haloes becoming increasingly aspherical 
and more prolate with increasing mass. 

The rounder haloes have less coherent rotation, with a 
median spin that is independent of mass (A mc( j(s > 0.9) « 
0.033). The most massive haloes have a median spin that 
is independent of sphericity (A me d(s) ~ 0.032). However, 
there is significantly less correlation between the nature of 
halo triaxiality (prolateness vs oblateness) and the spin pa- 
rameter. Although the haloes are far from being rotationally 
supported, there is a strong preference for the spin vector to 
be aligned parallel to the halo minor axis and to be perpen- 
dicular to the major axis. The tendency for the spin to be 
perpendicular to the intermediate axis is significantly weak- 
ened by the prevalence of prolate shapes for which there is a 
near degeneracy between the intermediate and minor axes. 

We find a clear signal that the spins and shapes of haloes 
are sensitive to the cosmological environment: more rapidly 
rotating haloes of a given mass are more strongly clustered. 
The strength of this effect increases with halo mass, ft is 
weak for subgalactic and galactic haloes, but can be larger 
than a factor of ~ 2 for galaxy cluster haloes. A similar 
effect is seen when examining halo shapes: more spherical 
haloes are more strongly clustered, with a greater signal at 



higher masses. Our result adds further evidence to the recent 
finding by Gao et al. (2005), also from analysis of the Mil- 
lennium Run, that the internal properties of haloes depend 
not only upon their mass but also upon the environment in 
which they form. 

The huge number of haloes in the ACDM Millennium 
simulation enables us to characterise the distribution of halo 
spins, and their relation to halo mass, shape and clustering, 
with unprecedented precision. However, we have also shown 
the significance of a careful halo definition. The properties of 
haloes defined and identified in different ways are noticeably 
different, and it is important to make the appropriate choice 
for a given application. For comparisons with real data, we 
recommend using the new class of "TREE" haloes which we 
have investigated in this work. 
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APPENDIX A: EXAMPLES OF HALOES 

The following plots give examples of haloes from our cata- 
logues, chosen as clear examples of various (usually undesir- 
able) groupfinder effects - they are not 'typical' haloes from 
their respective catalogues. We show them to give some vi- 
sual intuition as to the problems that can be encountered 
with different group-finding algorithms, as described in sec- 
tion 3. 

The first three show FOF haloes with unusual proper- 
ties in real and velocity space. Fig. Al shows a halo that is 
clearly made up of at least two objects joined via tenuous 
bridge in. We show a more massive halo that nevertheless 
consists of many linked objects in Fig. A2. Fig. A3 shows a 
very distorted object located near a much larger halo. 

The final two figures compare the results of the three 
different group-finding algorithms used. Figs A4 and A5 
compare haloes defined using the FOF, SO and TREE al- 
gorithms. In both cases, projections of the selected FOF 
halo (and its neighbours) are shown in the left-hand panels, 
and the corresponding SO/TREE halo and neighbours are 
in the right-hand panels. The selected haloes (green) have, 
again, been chosen to provide a striking illustration of the 
effects of different algorithms. The more 'normal' haloes in 
the background (blue) are less strongly affected by the choice 
of groupfinder. 
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Figure Al. An example of a FOF halo made of two objects linked by a tenuous bridge. The left-hand panels show projections of the 
halo (overplottcd in green) and its neighbours within 2 h~ Mpc (blue) in real space, in units of h _1 Mpc. The black contours show 
projected particle density, at 1, 10 and 100 particles per contouring bin. The potential- minimum centre of the main halo is marked with 
a red star, and the centre-of-mass is marked with a red cross. The right-hand panels show the particles of the selected halo in velocity 
space, in kms -1 . This halo has a mass of Mh = 11418m p 9.82 X 10 12 (i _1 M0, and a spin parameter of A = 1.5712. Its virial ratio, 
^ + 1 = —4.23 means it is excluded from the FOFclean catalogue. 
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Figure A2. As Fig. Al, but showing an example of a larger multi-object FOF halo. This halo has a mass of M h = 38 741m p 



3.33 X 10 h Mq, and a spin parameter of A = 0.3295. Its virial ratio, 44- 
catalogue. 



1 = —2.05 means it is excluded from the FOFclean 
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Figure A3. As Fig. Al, but showing an example of a small FOF halo with a large neighbour. The selected halo has a mass of 
Mh = 1967m p ft! 1.69 X 10 12 /i _1 Mq, and a spin parameter of A = 17.60. Its virial ratio, + 1 = —53.8 means it is excluded from the 
FOFclcan catalogue. 
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Figure A4. As Fig. Al, but showing a comparison of a massive FOF halo (left) and the corresponding SO halo (right), and their 
neighbours in each catalogue within 4fc _1 Mpc. A random sample of 1/4 of the particles are plotted; this doesn't effect the overall 
image of the halo. Projected logarithmic density-contours are plotted in black, spaced every factor of vlO between 1 and 10 3 ' 5 particles 
per contouring bin. The outermost contour (the halo boundary) is reproduced on the opposite panel in red. Note how the SO halo 
includes particles that were not part of the FOF group, and are much less dense than the halo proper. The FOF halo has a mass of 
M h = 744019m p ss 6.40 X 10 14 /i _1 M , and a spin of A = 0.05959. It has a virial ratio of ^ + 1 = -0.262, so it is included in the 
FOFclean catalogue. The SO halo has a mass of M h = 610 023m p pa 5.25 X 10 14 /i -1 Mq, and a spin of A = 0.04539. It has a virial ratio 



of 



1 = —0.332, slightly less relaxed than its FOF counterpart but still included in the SOclean catalogue. 
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Figure A5. As Fig. A4, but comparing a massive FOF halo (left) and the corresponding TREE haloes (right). Again, a random selection 
of 1/4 of the points are plotted, and the density contours are spaced every factor of vlO up to 10 2 ' 5 particles per contouring bin. Note 
how the FOF halo is highly extended, with the subhalo housing the potential minimum being linked to a large neighbouring halo. The 
TREE catalogue splits these into two objects. The FOF halo has a mass of M h = 126 033m p « 1.08 X lO 14 h _1 M , and a spin of 
A = 0.1953. It has a virial ratio of 2ZI -)- 1 = —0.111, so it is included in the FOFclean catalogue (despite its peculiar structure). The 
TREE halo has a much lower mass of M h = 40 719m p « 3.50 X 10 13 /i -1 M Q , and a spin of A = 0.05711. It is slightly less relaxed 
however, with a virial ratio of 33- + 1 = —0.159; it is included in the TREEclean catalogue. 



